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ABSTRACT 


A novel  mathematical  model  of  buoyant  convection  in  an  enclosure,  developed 
earlier,  is  solved  by  finite  difference  techniques  in  the  two-dimensional 
case.  This  model  has  been  developed  as  a principal  analytical  tool  for  the 
prediction  of  the  movement  of  smoke  and  hot  gases  in  fires.  Effects  of  large 
density  variations  caused  by  substantial  heating  are  retained  while  acoustic 
(high-frequency)  waves,  which  are  unimportant  to  buoyant  convection,  are 
analytically  filtered  out.  No  viscous  or  thermal  conduction  effects  are  in- 
cluded in  the  model.  These  two  characteristics  (filtering  and  no  dissipative 
effects)  distinguish  the  model  from  all  others  describing  buoyant  convection. 
The  mathematical  model  consists  of  a mixed  hyperbolic  and  elliptic  set  of  non- 
linear partial  differential  equations:  the  problem  is  a mixed  initial,  bound- 

ary value  one.  An  explicit  time-marching  algorithm,  second-order  accurate  in 
both  space  and  time,  is  used  to  solve  the  equations.  The  computational  proce- 
dure uses  a software  package  for  solving  a nonseparable  elliptic  equation 
developed  especially  for  this  problem.  The  finite  difference  solutions  have 
been  carefully  compared  with  analytical  solutions  obtained  in  special  cases  to 
determine  the  stability  and  accuracy  of  the  numerical  solutions.  The  computer 
model  has  been  used  to  compute  the  buoyant  convection  produced  in  an  enclosure 
by  a spatially  distributed  heat  source  simulating  a fire.  The  computed  results 
show  qualitative  agreement  with  experimentally  observed  buoyant  convection  in 
enclosure  fires. 
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I.  Introduction 


This  paper  presents  the  first  results  for  a finite-difference  integration  of 
an  approximate  set  of  equations  describing  buoyant  convection  in  an  enclosure. 
The  work  represents  a continuation  of  the  research  reported  in  Reference  1, 
where  the  set  of  approximate  equations  was  derived.  The  primary  application 
of  interest  to  the  authors  is  the  movement  of  smoke  and  hot  gases  caused  by  a 
fire  in  a room. 

The  research  presented  here  is  distinguished  from  previous  numerical  computa- 
tions of  buoyant  convection  in  three  respects.  First,  in  this  model  the  fluid 
is  taken  to  be  an  inviscid,  non-heat-conducting  perfect  gas,  and  the  spatial 
and  temporal  magnitude  and  variation  of  the  heat  source,  which  simulates  a 
fire  and  drives  the  flow,  are  taken  as  known.  These  approximations  are  justi- 
fied because  under  conditions  characteristic  of  even  a small  room  fire,  the 
Grashof  numbers  (representing  the  ratio  of  the  inertial  to  viscous  forces  for 
natural  convection)  are  large  enough  for  molecular  transport  phenomena  to  be 
important  only  in  wall  boundary  layers  and  in  the  highly  convoluted  flame 
sheets  which  constitute  the  region  of  intense  heat  addition.  The  study  of  the 
detailed  flame  structure  of  real  fires  is  an  extraordinarily  complicated  sub- 
ject in  its  own  right,  and  is  bypassed  here  by  specifying  the  heat  source. 

Wall  boundary  layers  represent  a local  refinement  to  be  considered  separately 
at  a later  date.  Batchelor  [2]  gives  a brief  but  relevant  discussion  of  the 
applicability  of  the  inviscid  equations  in  the  context  of  atmospheric  motions. 

It  should  be  noted  that  such  simplifications  do  not  preclude  a description  of 
turbulence;  but  no  turbulence  model  is  explicitly  included  in  this  study.  We 
note,  however,  that  any  turbulence  model  appended  to  the  present  equations 
must  be  of  the  "sub  grid"  variety,  since  no  spatial  or  temporal  averaging  is 
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Implied  in  the  equations  derived  in  Reference  1. 


The  omission  of  any  turbulence  model  is  based  on  the  observation,  quantified  by 
McCaffrey  [20],  that  most  of  the  energy  containing  fluctuations  in  buoyant 
plumes  induced  by  laboratory  scale  diffusion  flames  are  of  low  frequency  and 
large  spatial  extent.  Such  fluctuations,  with  frequencies  typically  in  the 
range  2 - 5 Hz  and  length  scales  comparable  to  the  local  plume  width  are 
directly  resolved  by  the  computational  procedure.  A knowledge  of  the  behavior 
of  higher  frequencies  and  smaller  length  scales  is  of  course  crucial  to  an 
understanding  of  combustion  related  phenomena.  However,  as  noted  above,  such 
questions  have  been  by-passed  in  the  present  study. 

Simple  models  of  smoke  and  hot  gas  transport  which  neglect  molecular  transport 
phenomena  have  been  reasonably  successful  in  predicting  global  properties  of 
flow  fields  [3]  . The  present  work  is  intended  as  a first  step  towards  more 
detailed  studies  along  these  lines. 

Second,  the  approximate  set  of  equations  integrated  in  this  paper  are  charac- 
terized by  the  fact  that  large  density  variations  due  to  temperature  changes 
are  admitted,  but  compressibility  effects  are  suppressed.  Such  a fluid  has 
been  termed  thermally  expandable  in  other  contexts  [4] . In  the  fire  setting, 
allowance  for  density  variations  due  to  temperature  increases  during  combus- 
tion is  essential.  It  is  common  for  temperature  in  a flame  to  exceed  1000°C, 
implying  that  the  density  decreases  locally  to  less  than  one-quarter  its  ambi- 
ent value  in  the  nearly  constant-pressure  process.  In  Reference  [1],  a set  of 
equations  of  motion  was  derived  formally  which  permit  description  of  large 
density  variations  in  a flow  while  ignoring  acoustic  oscillations  arising  be- 
cause of  the  elastic  properties  of  the  fluid.  Such  model  equations  include 
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the  important  features  of  buoyant  flows  without  requiring  excessive  computer 
time  necessary  to  determine  high-frequency  sound  waves  when  numerically  inte- 
grated. In  this  sense  the  equations  "filter  out"  the  sound  waves  while  de- 
scribing the  lower  frequency,  organized  motions  due  to  buoyant  effects  such  as 
internal  waves. 

Finally,  in  previous  efforts  to  compute  flow  fields,  produced  by  buoyancy  or 
by  any  other  mechanism,  there  are  few  reported  detailed  checks  on  the  quality 
of  the  numerical  solutions  to  the  finite  difference  equations.  By  contrast, 
in  References  8 and  22  detailed  comparisons  are  made  with  analytical/numerical 
solutions  obtained  to  the  general  difference  equations  in  simple,  special 
cases.  Through  these  comparisons  confidence  in  both  the  algorithm  and  its 
implementation  as  a computer  code  was  gained.  Such  tests  showed  that  the 
algorithm  is  stable  and  for  the  example  cases  performed,  the  error  made  in 
solving  the  difference  equations  at  any  time  step,  even  in  the  nonlinear  case, 
was  over  two  orders  of  magnitude  less  than  the  discretization  errors  made  in 
approximating  the  partial  differential  equations  by  finite  difference  equa- 
tions. For  buoyant  flows  of  the  type  considered  in  this  paper,  it  is  espe- 
cially important  to  have  confidence  in  both  the  stability  and  the  accurancy 
of  the  algorithm  so  that  the  real  physical  instability  represented  by  the 
fluid  motion  can  be  distinguished  from  any  computational  instabilities. 

In  Section  II  the  equations  derived  in  Reference  1 are  recast  into  the  form  in 
which  they  are  solved,  and  the  finite  difference  approximations  to  the  equa- 
tions are  presented.  In  Section  III  solutions  determined  by  this  model,  for 
the  buoyant  flow  produced  by  a heat  source  in  two  rectangular  enclosures,  are 
presented  and  discussed. 
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II.  Formulation 


A.  Continuous  Problem 

In  Reference  1 the  authors  derived  a set  of  nonlinear  equations  describing 
nondissipative , buoyancy  driven  flows  of  a perfect  gas.  The  flows  were  assumed 
to  be  generated  by  a localized  heat  source  in  which  the  heat  is  added  slow- 
ly so  that  the  time  scale  associated  with  the  heat-source  grov;th  and  resultant 
fluid  motion  is  long  compared  with  the  transit  time  of  an  acoustic  signal 
across  the  spatial  extent  of  the  source.  Flows  Induced  by  a room  fire  gener- 
ally satisfy  this  assumption.  Properties  of  the  equations  were  discussed  in 
Reference  1.  In  this  section  these  equations  are  rewritten  in  a form  appro- 
priate for  numerical  integration  by  finite  difference  techniques,  and  the 
boundary  conditions  for  the  equations  are  presented. 

As  in  Reference  1 we  consider  an  inviscld,  non-heat-conducting  perfect  gas. 

The  magnitude  and  the  spatial  variation  of  the  heat  source  (representing  the 
exothermic  reaction  in  a fire)  are  taken  as  known;  justification  for  such  a 
model  is  given  in  Reference  1.  The  fluid  and  the  fire  source  are  assumed  con- 
fined in  a closed  rectangular  room  with  the  center  of  the  source  along  the 
floor.  In  contrast  to  Reference  1,  we  consider  only  a completely  enclosed 
room  (no  leaks),  and  when  difference  equations  are  Introduced,  we  confine 
attention  to  the  two  dimensional  evolution  of  the  flow. 

The  continuity,  momentum,  energy,  and  state  equations  are  given  respectively 
in  Reference  1 as: 


A 


8p  3 

— + (puj[)  = 0 

9t  3xi 


3ujl  9u^  3(p-pQ(t)) 

+ U4  ) + pk^g  = 0 

3t'’3xj/  3xj^  (1) 


/ 3T  3T  \ 

( — + ) Q(xi,t) 

^ 3t  ^ 3xj/  dt 


Po(t)  = pRT 

th 

Here  p is  density,  Uj^  the  velocity  in  the  i coordinate  direction 
(i  = 1,  2,  3),  p is  the  pressure  excess  above  the  mean  pressure  Po(t)  in  the 
room,  T the  temperature,  Cp  the  constant-pressure  specific  heat,  R the  gas 
constant,  k^g  is  the  gravitational  acceleration  (of  magnitude  g)  and  Q(xj^,t) 
the  specified  volumetric  heat  source*  The  spatially  uniform  mean  pressure 
PqCp)  depends  only  upon  time  and  increases  because  of  the  heating  within  the 
the  room.  It  is  determined  in  a completely  enclosed  room  by  the  equation 

^Po  y-1 

= / 0(xi,t)dV  (2) 

dt  V V 

where  y is  the  ratio  of  specific  heats,  V is  the  volume  of  the  room  and  the 
integration  is  performed  over  this  entire  volume.  Equation  (2)  is  a thermo- 
dynamic statement  that  the  mean  pressure  rise  as  a function  of  time  is  deter- 
mined by  the  total  heat  added  to  the  room.  (Heat  can  only  be  added  or  removed 
volumetrically  and  not  through  the  walls  because  thermal  conduction  and  radi- 
ative transport  have  been  ignored  in  this  model.)  It  will  also  turn  out  to  be 
a mathematical  consistency  condition  required  if  a solution  for  the  pressure 
field  is  to  exist. 
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Equations  (1)  and  (2),  are  the  approximate  set  of  nonlinear  equations  solved 
by  finite  difference  techniques  in  two  spatial  dimensions  in  this  paper.  The 
equations  admit  buoyant  or  internal-wave  motions  while  "filtering  out"  high- 
frequency,  acoustic  waves.  They  reduce  to  the  Boussinesq  equations  when  heat- 
ing is  mild,  total  density  variations  are  small,  and  variations  in  the  mean 
background  pressure  can  be  neglected  (as  would  be  the  case  if  the  room  con- 
sidered here  were  open  or  if  the  mean  pressure  variation  were  comparable  to 
the  spatial  pressure  perturbation) . To  recast  the  equations  into  a form  more 
suitable  for  numerical  computation,  we  take  the  substantial  derivative  of  the 
equation  of  state  and  use  this  with  the  energy  equation  to  eliminate  the 
temperature.  The  resulting  equation  describes  the  evolution  of  the  density 
under  heating 

3p  3p  3^i 

— + U£  = - p = - pD(x-;,t)  (3) 

3t  3xj[^  3x^ 


where 


D(xj ,t) 


1 

YPo(t) 


[ 


(Y-l)Q(xj  ,t)  - 


(4) 


Equation  (3)  and  the  continuity  equation  identify  D(x^,t)  as  the  divergence 


3uj[ 

= D(x4 ,t)  (5) 

3xj^  ^ 


Finally,  as  in  Reference  1,  the  equation  for  the  spatially  variable  portion  of 
the  pressure  is  obtained  by  dividing  the  momentum  equations  by  density  and 
taking  the  divergence  of  these  equations.  The  resulting  equation  is 

3 

3x£ 


3 / 1 3p 

Sxj’  'p  3x£ 


3ui  , 3D(x-j  ,t) 

+ — __ 


(6) 


3x4 


3t 
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Equation  (6)  is  the  generalization  for  a "thermally  expandable"  fluid,  which 
we  consider  here,  of  the  well-known  incompressibility  condition  in  Boussinesq 
fluids.  When  the  density  is  constant,  then  D(x,t)  = 0,  and  Equation  (6)  re- 
duces to  Poisson's  equation.  The  boundary  conditions  on  these  equations  are 
that  velocity  normal  to  any  (impermeable)  wall  vanish. 

u^nji^  = 0 (7) 

where  n;j^  are  the  normal  components  of  a vector  describing  the  boundary  walls. 
From  Equations  (1)  and  these  conditions,  the  appropriate  boundary  conditions 
on  the  pressure  equation  are  obtained 

3p 

n^  = pg  n^kj^  (8) 

8ni 

An  important  observation  must  be  made  concerning  Equation  (6)  and  Neumann 
boundary  conditions  (8).  When  Equation  (6)  is  integrated  over  the  total  vol- 
ume of  the  room,  both  the  left  hand  side  of  the  equation  and  the  first  terra  on 
the  right  are  divergence  forms  and  can  be  converted  into  integrals  over  the 

boundaries  of  the  room.  Application  of  boundary  conditions  (7)  and  (8)  show 

3D 

that  each  of  the  terms  is  zero;  therefore,  the  integral  over  the  volume  of  — 

3t 

must  also  be  zero.  The  requirement  that  the  integral  of  Equation  (4)  for  D 
over  the  volume  be  zero  produces  Equation  (2),  the  condition  for  the  spatially 
uniform  background  pressure.  Therefore,  the  elliptic  equation.  Equation  (6), 
for  the  pressure,  with  the  Neumann  boundary  conditions,  is  seen  to  produce  a 
condition  which  must  be  satisfied  for  a solution  to  the  equation  to  exist. 

This  condition,  Equation  (2),  determines  the  time  evolution  of  the  spatially 
uniform  background  pressure  and  demonstrates  that  the  total  pressure  can  be 
consistently  separated  into  a spatially  uniform  background  pressure  and  an 
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inhomogeneous  time  dependent  over-pressure.  In  the  next  section  describing 
the  difference  equations,  exactly  analogous  considerations  are  found  to  apply 
to  the  linear  algebraic  equations  approximating  Equation  (6)  and  boundary 
conditions  (8)  . 

For  selecting  a difference  scheme,  the  second  of  Equations  (1),  the  momentum 
equations,  are  rewritten  in  vector  invariant  form,  noting  that  uj[  are  corapo- 
nents  of  the  velocity  vector  field  u 

->• 

1 -»• 

— + 1/2  V (q2)  - uxto  = --  V p+kg  (9) 

9t  p 

where  q^  = u*u  and  w = V x u is  the  vorticity.  The  curl  of  Equation  (9) 
yields  the  vorticity-transport  equation 


^ /lx 

— - Vx(uxo))=-V  (-)  xVp 
3t  ' P ^ 


(10) 


Since  u is  a vector  field,  it  is  necessary  to  calculate  correctly  both  its 
divergence  and  its  curl.  Equation  (5)  specifies  the  divergence  of  u,  and  the 
equation  for  the  pressure.  Equation  (6),  assures  that  the  divergence  is 
properly  determined  at  each  time.  Equation  (10)  is  the  equation  for  the  evo- 
lution of  the  curl  of  the  velocity,  the  vorticity.  The  difference  scheme 
selected  must  assure  that  Equation  (10)  is  satisfied  in  difference  form. 

The  complete  set  of  recast  nonlinear  equations  are  gathered  and  rewritten 
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below 


9p  9p 

— + U4  = - pD(x-  t) 

9t  9xi 


where 


dpo 

dt 


Hi 

V 


/ Q(xi^t)  dV 
V 


D(xi^t) 


1 

YPo(t) 


dPo 

(y-l)Q(xi^t) 


is  the  permutation  tensor  and 


9% 


o>i  = Ei-jii  ate  the  components  of  the  vorticity  vector. 

9xj 


Boundary  conditions  are 


Uini  = 0 


ni  — - = p g ki  ni 
9xj^ 


(lla) 

(llb) 

(llc) 

(lid) 

(lie) 


(12a) 

(12b) 
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For  numerical  integration,  Equations  (11)  and  boundary  conditions  (12)  are  al- 
tered in  two  additional  ways.  First,  we  calculate  the  density  and  pressure 
differences  from  their  initial  values,  which  may  be  functions  of  the  vertical 
coordinate.  This  is  done  to  eliminate  accuracy  problems,  since  the  thermally 
induced  density  and  pressure  differences  can  be  very  small  during  the  early 
portion  of  the  heating  process.  Second,  the  equations  are  made  dimensionless. 
The  nondimenslonalization  is  done  both  for  convenience  and  to  ensure  proper 
scaling.  All  dependent  quantities  are  made  to  be  of  order  unity  in  magnitude 
for  purposes  of  the  computation.  The  remaining  dimensionless  parameters 
characterize  the  strength  and  location  of  the  heat  source  as  well  as  the  room 
geometry. 

In  Figure  1 a schematic  diagram  of  a fire  evolving  a room  and  a set  of  coordi- 
nate axes  are  shown.  It  is  assumed  that  initially  the  enclosure  is  filled 
with  quiescent,  stratified  fluid  of  density  Po(y)»  where  we  denote  xj  = x, 

X2  = y and  X3  = z.  We  define  a density  difference  from  ambient  and  a pressure 
difference  as  follows: 

p (x,y,z,t)  = p(x,y,z,t)  - Po(y)  (13a) 

y 

p (x,y,z,t)  = p(x,y,z,t)  - Po(t)  + g / Po^y*)  ‘^y'  (13b) 

o 

These  differences  p and  p need  not  be  small  compared  with  Po(y)  ^ad  Po(t) 
respectively.  Then  Equations  (lla)-(llc)  become 

3p  3p  dpQ 

— + Uj[ + V = - (Po(y)  + p)  D(x:j  t)  (14a) 

3t  8xi  dy  ’ 
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= " (l/(Po+p)) 


(lAb) 


3ui  3 

+1/2  (u.u.)  - 

at  axi 


9p 

3X;{^ 


- k^gp 


(” 
V 8x£ 


~ 9p  \ 

(l/(Po+p)) ) 

ax^/ 


a ~ ~ 3D  3^  . 

— (gp/(Po+P))  + — + 1/2  r (q^) 

8y  . at  3x^2 


3 

8xj[ 


where  v = U2  and 

and  the  boundary  condition  (12b)  becomes 


(14c) 


3p 

Hi  = p gki  n^  (15) 

8x£ 


Finally,  we  form  dimensionless  equations  using  the  density  Pqq  = Pq(0),  the 
height  of  the  room  H as  the  length  scale  and  the  free  fall  time  (H/g)^/^  as 
the  time  scale.  Then,  denoting  dimensionless  quantities  with  a hat 

A.  A »N* 

P ~ P/PoO  » P “ p/poo^n  » pQ  ~ Po'^PoO 

(16) 

A A 

u^  = u^/ZgH  , Xji^  = Xj^/H  , t = t//  (H/g) 


Equations  (13)  remain  exactly  the  same  in  dimensionless  form  with  g set  equal 
to  one.  Subsequently,  in  this  paper  all  quantities  will  be  understood  to  be 
dimensionless,  and  the  hat  notation  will  be  dropped.  For  the  dimensionless 
coordinates,  we  note  that  0 ^ x l/AR  , 0 ^ y ^ 1 and  0 ^ z _<  1/BR  where 
AR  = H/L  and  BR  = H/W,  Also,  in  the  remainder  of  this  paper,  the  problem  will 
be  specialized  to  two  spatial  dimensions  so  that  all  quantities  will  be  as- 
sumed independent  of  z. 
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This  non-dimenslonallzation,  while  simplifying  the  form  of  the  equations,  does 
not  relate  the  scale  of  the  induced  motion  to  that  of  the  source.  An  alter- 
native scheme  which  does  have  this  feature  is  derived  in  the  Appendix.  The 
resulting  equations  are  somewhat  less  convenient  for  numerical  computation, 
except  in  the  Boussinesq  limit  when  density  differences  are  small.  Hence, 
the  variables  as  defined  by  Equation  (16)  will  be  used  throughout  the  remain- 
der of  the  paper. 

B.  Finite  Difference  Equations 

1.  Basis  for  the  Selection  of  the  Difference  Equations 
In  this  section  the  finite  difference  equations  and  the  boundary  relations  for 
the  solution-algorithm  are  presented.  This  algorithm  does  not  represent  a new 
numerical  scheme;  rather  it  represents  an  adaptation  of  several  well-known  and 
proven  ideas  to  the  equations  of  interest  here.  Since  these  equations  are 
different  than  those  generally  examined  in  computational  fluid  dynamics,  there 
are  interesting  features  which  vrill  be  noted  as  the  difference  equations  are 
presented . 

It  is  important  first,  however,  to  state  the  requirements  which  we  used  in 
selecting  the  scheme.  These  requirements  do  not  necessarily  specify  a unique 
scheme,  but  restrict  greatly  the  selection.  As  noted  in  the  previous  subsec- 
tion, the  momentum  conservation  equations  provide  relations  for  determining 
the  velocity  field,  which  is  a vector  field,  as  a function  of  time.  Forming 
the  difference  equations  from  the  vector  invariant  form  of  the  momentum  equa- 
tions is  one  criterion.  This  choice  is  made  to  assure  that,  when  the  discrete 
analogues  of  the  divergence  and  curl  are  applied  to  the  difference  equations. 
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suitable  discretized  equations  for  the  pressure  and  the  vorticity  transport 
are  obtained.  Thus  we  can  be  certain  that  the  divergence  and  the  curl  are 
calculated  correctly  in  discrete  form,  the  divergence  producing  the  equation 
for  the  pressure  and  the  curl  producing  an  evolution  equation  for  vorticity. 

A second  criterion  imposed  is  that  the  difference  equations  provide  second- 
order-accurate  approximations  to  the  partial  differential  equations.  This 
condition  is  imposed  because  second  order  accuracy  is  necessary  to  obtain 
reasonable  spatial  and  temporal  resolution  for  the  hundreds  of  time  steps  re- 
quired to  calculate  the  complete  evolution  of  the  flow  field  in  a room  fire. 

A well-understood  buoyant  flow  field  is  that  produced  by  internal  waves  in  an 
ambient  stratified  environment.  Internal  waves  have  been  analyzed  and  calcu- 
lated for  a long  time  [6,7];  they  are  determined  as  solutions  to  the  linear- 
ized partial  differential  equations  of  buoyant  flow.  In  addition,  internal 
waves  can  be  expected  to  arise  naturally  in  a room  fire  setting  when  the 
driving  fire  has  heated  and  stratified  the  room  air;  after  the  fire  has 
extinguished  itself  because  of  vitiated  air,  for  example,  pure  internal  wave 
modes  will  exist.  An  additional  criterion  which  was  imposed  on  the  selection 
of  the  difference  scheme  was  that  it  accurately  reproduce  internal-wave  modes 
in  an  enclosure:  if  the  difference  scheme  does  not  reproduce  this  rather 

simple  linear  flow  field,  it  is  unlikely  to  reproduce  more  complicated  flows. 
In  Reference  8 the  authors  examined  internal  waves  and  some  second  order 
linear  difference  equations  that  reproduce  these  waves.  This  analysis 
determined  the  difference  scheme  for  all  but  the  non-linear  convective  terms 
in  the  momentum  equations. 
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In  an  important  paper  for  numerical  weather  forecasting  [9],  Arakawa  discusses 
design  of  computational  schemes  for  long-term  numerical  integration  of  the 
equations  of  fluid  motion  in  the  case  of  two-dimensional  incompressible  flow. 
The  major  thrust  of  his  paper  is  a derivation  of  second-order  accurate  spatial 
difference  schemes  which  eliminate  the  nonlinear  computational  instability 
first  noted  by  Phillips  [10].  Arakawa  emphasizes  that  for  two-dimensional, 
incompressible  flow,  the  discrete  approximations  to  quadratic  forms  of  de- 
pendent variables,  such  as  the  velocity  squared  or  the  vorticity  squared  (or 
both),  must  be  conserved  when  the  continuous  variables  are,  and  he  uses  this 
constraint  to  select  three  difference  approximations  which  are  acceptable. 
Another  criterion  imposed  in  the  selection  of  our  difference  scheme  is  that 
it  reduce  to  one  of  Arakawa* s acceptable,  or  stable  schemes  with  respect  to 
the  nonlinear  computational  instability  when  the  flow  is  incompressible.  The 
scheme  chosen  is  denoted  J3  by  Arakawa  [11]  and  conserves  energy  in  the  in- 
compressible case.  It  can  be  obtained  by  differencing  the  vector  invariant 
form  of  the  momentum  equations. 

Finally,  we  v;anted  the  difference  scheme  to  be  easily  generalized  to  three- 
dimensional  flow  configurations.  For  a model  of  buoyant  flow  driven  by  a fire 
to  be  successful,  it  must  be  able  to  calculate  three-dimensional  flows.  The 
scheme  selected  and  presented  below  satisfies  all  of  the  criteria  stated 
above . 


2.  The  Equations 

In  Figure  2a,  the  two-dimensional  rectangular  enclosure  in  dimensionless 
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variables  is  shown  together  with  a schematic  representation  of  the  spatial 

grids  used  for  the  finite  difference  scheme.  The  grid  formed  from  solid  lines 

represents  the  basic  mesh  into  which  the  enclosure  is  divided:  in  general 

there  are  I mesh  cells  in  the  x-direction  and  J mesh  cells  in  the  y-direction. 

Upon  this  basic  mesh,  the  two  components  of  the  vector  velocity  (u,  v)  and 

8 V 3u 

single  surviving  component  of  the  vector  vorticity  w = — - — are  defined. 

3x  8y 

The  second  grid,  formed  by  joining  the  center  points  of  the  basic  grid  cells 
and  denoted  by  dashed  lines,  is  that  upon  which  scalar  quantities  such  as  den- 
sity p and  pressure  p are  defined.  In  Figure  2a  the  densities  in  the  left- 
hand  column  of  cells  and  in  the  bottom  row  of  cells  are  shown  to  indicate  how 
they  are  enumerated  for  the  numerical  computation. 

In  Figure  2b  a typical  mesh  cell  is  shown.  Illustrating  where  all  of  the 
dependent  variables  in  the  finite  difference  scheme  are  defined  relative  to 
the  cell. 

The  following  discretely  evaluated  functions  will  denote  approximations  to  the 
corresponding  solutions  to  Equations  (9)  and  (10); 

n 

u = u(i6x, ( j-l/2)6y ,n5t) 

ij 

V = ((l-l/2)6x,  j6y,  n6t) 

Ij 

p = p((i-l/2)6x,  (j-l/2)6y,  n6t) 

ij 
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p = p((i-l/2)6x,  (j-l/2)6y,  n6t) 

ij 


D = D((i-l/2)6x,  (j-l/2)6y,  nfit)  , 

ij 


n 

w = w(i5x,  i6y,  nSt) 

ij 


where  6x  = 1/(I»AR)  and  6y  = 1/J  are  the  mesh  cell  sizes  in  the  x-  and  y- 
directions  respectively  and  where  6t  is  the  time-step  size.  Such  a staggered 
grid  is  commonly  used  for  multidimensional  finite  difference  integrations 
[12]  . 


With  this  notation,  the  following  set  of  finite  difference  equations  was  used 
to  approximate  Equations  (11)  and  boundary  conditions  (12): 


For  Equation  (14a),  l<i<I,  l<j<J  and  n > 2, 


-n+1 

) 

ij 


1 


1+(1/2)D  6t 

ij 


!~n-l  / \ 

p (l-(l/2)D  6t) 

ij  V ij  / 


~26t 


n n n ■ 

F + F + (1/2)D  Po(j) 
pxij  pyij  ij 


where 


(18) 


n 

p = p - Po(j)  = the  density  difference  from  the  initial  density, 

ij  ij 


Po(j)  = exp[-( j-l/2)6y/Yg]  = the  prescribed  ambient  initial  stratification, 

(19) 


Yg  = the  stratification  length  scale. 
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n n 

The  flux  terras  F and  F for  l<i<I,  l<j<J  are  given  by 

pxij  P7ij 


n 

F 

PXij 


n 


~n  ~n  \ / n “ 

i,j+l  ij  i.j-1 


pyij 


6y 


~n  ~n  \ / n n 

Po(3+l)  - + P “P  \ / ^ 

i,j+l  ij  l>j-l 

26y 


(20b) 


Equation  (18)  employs  a modification  of  the  second-order  accurate  central  dif- 


ference (leap  frog)  temporal  discretization.  The  modification  eliminates  an 


instability  that  would  arise  if  the  leap  frog  scheme  had  been  applied.  It 


affects  the  undifferentiated  term  pD(x,y,t)  in  equation  (14)  that  is  well 


known  to  lead  to  a computational  instability  for  ordinary  differential  equa- 


tions when  leap  frog  differencing  is  used  [13], 


For  Equations  ( 14b) 


n+ 1 n- 1 

u = u - 26t 

ij  ij 


P - P 

i+l,j  ij 


(21a) 


.n 


2p„(  j ) + P + P 

i+l,j  ij 
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and  for  KKI-l,  Kj<J  , 


n+1 

n-1  ' 

n 

2 

V 

= V - 26t  < 

F 

+ — 

ij 

ij  I 

•H 

6y 

(21b) 


for  Ki<I,  Kj<J-l  . 


n n 

The  fluxes  F and  F are  defined  as  follows; 
xij  yij 


for  KKI-I,  Kj<J 


n 1 


xlj  26x 
and  for  Kj<J-l,  Ki<I 


n \2  2 


1 / n n n n 

V 03  + V 03 

2 \ 03 . , ij  w . . , i , 

ij  i,J-l  ^ 


(22a) 


n 1 

yij  26y 


n 


i.i+1 


n n n 

U 03  + u 

03 


n 


03 


ij  ij  ‘‘^i-l.j 


(22b) 


and  where 


n n 


n 


n 


ij 


V. , , . - V. . U.  . , 1 - u. . 

i+l,j  ij  i,j+l  ij 


5x 


6y 


n In  n n ln  n 

V =-(v  +v  ) ,u  =-(u  +u) 

03  2 ij  i+l,i  Wj.  2 i»j+l  ij 


ij 


ij 


(22c) 


n 2 

(q  ) = 

ij 


n n 

u. . + u.  , . 

ij  i-l,J 


n n 

V. . + V. 


ij  i.j-1 
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Note  that  boundary  conditions  (12a)  on  the  normal  velocities  imply  that 

UQ^j  = =0  for  l<jO  and  = 0 for  Ki^I»  These  boundary  condi- 

n n 

tions  are  applied  formally  in  the  expressions  for  the  fluxes  F , F , 

pxij  pyij 

n n 

F and  F in  mesh  cells  adjacent  to  boundaries, 
xij  yij 


The  finite  difference  analog  to  Eq,  (14C)  is,  for  Ki_<I  and  KjO 


2 

p 

i+l,j 

~n 
• P 

ij 

^n 

P 

ij 

~n 

P 

i-l,j 

6x^ 

^n 

^n 

~n  ^n 

2Po(j)  + P 

+ P 

2po(j)  + P 

+ P 

^ i+1 , j ij 

ij  i-l»j-' 

2 


+ 


6y2 


~n  ^n 

P “ P 

i,j+l  ij 


~n  ^n 

Po(j+l)  + PoO)  + Pi  j+i  + Pij 


P - P 

ij  i,j-l 


~n 

Po^j)  + Po(j-l)  + Pij  + Pi,j-1  - 


n+1  n-1  n n n n 

D-n  F -F  F F 

ij  ^°ij  x,i-l,j  x,ij  yi,j-l  “ yij 
+ + 

26t  6x  5y 


1 


<Sy 


^n  ^n 

P + P 

i.j+1  ij 


~n  ^n 

Po(j+l)  + Po^j)  + P + P 

i,j+l  ij 


^n  ^n 
P + P 

ij  i»j-l 


~n  ^n 

PoO)  + PoO-1)  + P + P 

ij  i,j-l'^ 


(24) 
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This  difference  equation  for  the  pressure  arises  formally  by  applying  the 
finite  difference  analogue  of  the  divergence  operator  to  Equation  (21)  and 
noting  that  the  finite  difference  divergence  of  the  velocity  field  satisfies 
the  equation 

n n n n 

^Ij  " ^i-l,j  _ ^ij  " ^i,j-l 

-J-  ■ ■ ' ■ 

6x  6y 

(This  equation  is  the  difference  approximation  to  Equation  (5).) 

The  boundary  conditions  (15)  in  discrete  form  become 


= D (25) 

ij 


~n  ^n 


P = P 
0,j  l,j 


for  KjO 


~n  ^n 
P = P 

1+1,  j 


(26a) 


~n 

P = 5y 
i.O 


^n  ^n 
P + P 
. i,l  i,0 


/2 


for  Ki<I 


~n 

/ ~n 

~n  \ 

P 

- p = - 5y 

P + P 1 /2 

1,J+1 

1,J 

V i,J+l 

i,J  / 

(26b) 


~n  ~n 

Although  p and  p appear  in  the  boundary  conditions,  Equation  26(b), 

i,0  i,J+l 


they  also  appear  in  Equation  (24)  with  j = 1,J  in  the  same  combination.  As  a 


~n  ~n 

result,  p and  p never  need  to  be  specified  to  obtain  a solution  to 

1,0  i,J+l 
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Equation  (24)  and  the  boundary  conditons  (26).  Equations  (24)  together  with 

boundary  conditions  (26)  constitute  a singular  linear  algebraic  system  of 

equations.  When  Equations  (24),  with  boundary  conditions  (26)  incorporated, 

are  summed,  the  left  hand  side  sums  to  zero,  demonstrating  that  all  of  the 

equations  are  not  linearly  independent.  Also,  the  last  three  terms  on  the 

right  hand  side  sum  to  zero,  producing  the  requirement  that  the  double  sum 
n+1  n-1 

Dij  - Dij 

over  must  vanish. 

26t 

The  vanishing  of  the  sura  of  the  left  hand  sides  of  Eqs.  (24)  and  (26)  is  the 
discrete  analog  of  the  fact  that  the  left  hand  side  of  the  integral  of  Eq.  (6) 
over  the  room  volume  vanishes  when  the  boundary  conditions  u • n = 0 are 
applied.  The  vanishing  of  the  sum  of  the  right  hand  sides  of  Eqs.  (24)  and  (26) 
is  the  corresponding  discrete  analog  for  the  requirement  on  the  integral  of  the 
right  hand  side  of  Equation  (6).  This  requirement,  that  '['l  must  vanish  at 
each  time  level,  is  the  discrete  analog  of  Eq.  (lid). 

n 

Examination  of  Equation  (28)  for  shows  that  it  has  been  chosen  so  that  its 

double  sum  over  all  mesh  points  vanishes,  and  that  the  condition  which  must  be 

satisfied  to  allow  this  choice  produces  Equation  (29)  for  the  mean  pressure. 

Therefore,  the  singular  linear  algebraic  system  is  seen  to  be  consistent  and 

thus  to  allow  a solution.  The  solution  is  made  unique  by  specifying  that  the 

~n 

double  sum  over  all  mesh  points  of  pj^^j  is  zero.  This  is  tantamount  to 

~n 

specifying  that  Pq  is  literally  the  mean  pressure  in  the  room,  with  p^j  the 

perturbation  about  the  mean.  Details  of  the  algorithm  used  to  solve  Equations 

~n 

(24)  and  (26)  for  pj^j  are  presented  in  Reference  17. 
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The  heat  source  has  been  chosen  to  have  the  form 


n n 

Q = 0 f 

ij  ij 


Q = {-)  X exp  [-3(xi“Xc)^  ~ Xy^] 

ij  ^ 


(27a) 


(27b) 


xi  = (i-l/2)6x  , Yj  = (j-l/2)6y 


(27c) 


n n 

f = Qo  tanh  At 


o 

t = 0 


n 

t 


n-1 


I 

n * =0 


Hence,  the  discrete  divergence  of  the  velocity  field  becomes 

^1  " n 

D = [(Y-l)Qii“K]f 

ij  YPo” 


(27d) 


(27e) 


(28a) 


where 


K = 


Y-1 

IJ 


I 


I 

i=l 


J - 


I Qij 

j=i 


(28b) 


and  the  mean  background  pressure  is  found  from  the  difference  equation 


n+ 1 n-1  n n 

p = p + Kf  25t 
o o 


(29) 


0 1 o 

with  p = p =1  since  f = 0 . 
o o 
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Since  difference  equations  (18)  and  (21)  are  three-level,  second-order  schemes 
(leap-frog)  in  time,  a starting  procedure  is  needed.  The  following  first- 
order,  explicit  scheme  was  used  to  start  the  computation  and  to  restart  when  a 
change  in  the  time  step  has  been  made: 


n+1  n 

1 

L = u - 6 1 

{ ^ 

ij  ij 

1 xij 

i+l,j 


-n 


2Po(j  ) + P + 

i+l,j 


(30a) 


n+1 

n 

n 

2 

= V 

- 6t 

F 

+ — 

ij 

ij 

yij 

6y 

^n  ^n  /^n  ^n  \ 6y 

P -P  + P +P  — 

iJ+1  ij  \ i,j+l  ij/  2 


~n  ^n 

Po(j+0  + Po^j)  P P 

i,j+l  ij 


(30b) 


~n+l  ~n 
P = P 

ij  ij 


~n 

Po(j)  + P 

ij 


n 

D 

ij 


n 

+ F 

Pxij 


n 

+ F 


Pyij 


(31a) 


n+1  n n 

p = p + Kf  6t 
o o 


(31b) 


n+1  ~n+l  ~n 

When  starting.  Equations  (31)  are  used  to  obtain  p and  p . Then  p 

° ij  ij 


n+1  n-1 

Dij  - 

is  obtained  from  Equation  (24)  with  replaced  by 

25t 


n+1  n 
i*ij  “ i>ij 

5t 
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~n  n+1 

With  this  solution  for  p , Equations  (30)  are  used  to  obtain  u and 

n+1 

V . The  starting  procedure  is  completed  (and  two  levels  of  all  dependent 

ij 

variables  have  been  obtained)  by  solving  Equation  (24)  with  n replaced  by 
n+1  throughout.  Subsequent  time  steps  are  taken  in  a straightforward  fashion 
with  the  density  and  velocity  components  being  advanced  through  Equations  (18) 
and  (21),  and  the  pressure  being  updated  through  Equation  (24). 

The  linear  stability  of  the  algorithm  is  the  only  other  consideration  for 
discussion.  A linear  stability  analysis  of  Equation  (18)  for  the  density 
shows  that  the  time  step  6t  must  satisfy  the  following  condition  for 
stability 


where 


n 

2 

6t  £ 

max 

+ 

KKI 

(V  ij/ 

6x 

Kj<J 

n 

n 

\ 

U 

= (1/2) 

(xx  + 

u 

) 

ij 

V ij 

i-1 

»j  / 

n 

n 

V 

= (1/2) 

(v  + 

V 

) 

ij 

^ ij 

i.j 

-1^ 

6y 


(32) 


When  the  stability  condition,  Equation  (32),  is  not  satisfied  by  a time  step, 
n 

the  time  step  6t  is  halved.  Then  the  time-marching  algorithm  is  restarted 


using  the  last  time-level  values  as  Initial  conditions.  A first-order  time 


step  is  taken  and  then  leap  frog  is  resumed. 
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When  the  finite  difference  analogue  of  the  curl  (see  Equations  (22))  is 
applied  to  the  difference  equations  for  the  velocity  components,  Equations 
(21),  a discretized  form  of  the  equation  for  the  vorticity  transport.  Equation 
(10),  is  formed.  A linear  stability  analysis  of  the  difference  equation 
yields  exactly  the  same  form  for  the  stability  criterion  as  that  found  above 
for  the  density  equation.  Reference  to  Figure  2b  shows  that  the  density  and 
vorticity  are  evaluated  at  different  points  in  the  mesh,  however,  and  there- 
fore, the  divergence  D and  the  velocity  components  U and  V are  to  be  evaluated 
at  different  points  than  those  used  in  Equation  (32),  To  account  for  the  dif- 
ference in  the  stability  criterion  implied  by  the  different  mesh  location 
points,  in  all  calculations  performed  using  the  algorithm  described  above,  the 
time  step  was  chosen  to  be  less  than  or  equal  to  0.8  the  maximum  value  found  for 
the  right  hand  side  of  Equation  (32). 

III.  Example  Calculations 

The  algorithm,  described  in  Section  II  and  tested  as  discussed  in  Reference  22, 
has  been  used  to  compute  solutions  to  the  buoyant-flow  equations.  In  this 
section  results  of  two  computations  are  presented  and  discussed.  One  calcula- 
tion is  for  the  flow  generated  by  a heat  source  centered  along  the  floor  in  a 
square  enclosure,  and  the  other  is  for  the  flow  generated  by  a heat  source  with 
maximum  along  the  floor,  but  off-center  in  a rectangular  enclosure,  twice  as 
long  as  it  is  high.  Other  calculations,  exploring  parametrically  the  features 
of  this  model  and  examining  in  detail  the  numerical  results  by  analyzing  the 
flow  data  computed  by  the  model,  will  be  reported  in  a companion  paper. 


25 


A.  Square  Enclosure  with  Centered  Heat  Source 
In  Figure  3 contours  of  constant  temperature  (isotherms)  are  shown  at  dimen- 
sionless time  2.0  for  a volumetric  heat  source  centered  along  the  floor  in  a 
square  room.  The  rate  of  heat  added  per  unit  volume  is  largest  along  the  floor 
at  the  center  of  the  room  and  decreases  in  a Gaussian  fashion  with  horizontal 
distance  from  the  center  and  exponentially  with  height  above  the  floor:  the 
dependence  of  the  heat  source  upon  position  in  the  room  is  given  by  Equation 
(27)  with  = 0.5.  The  heat  source  is  "turned  on"  slowly  according  to  Equa- 
tions (27a)  and  (27d)  asymptoting  to  full  strength  around  t = 10.0.  At  this 
early  time  the  problem  is  still  linear;  the  flow  velocities  are  sufficiently 
small  that  convection  is  unimportant,  and  the  temperature  increase  in  the  fluid 
is  directly  proportional  to  the  volumetric  rate  of  heat  added.  Therefore,  the 
isotherms  are  also  contours  along  which  the  volumetric  heat-addition  rate  is 
constant.  (These  contours  can  be  seen  to  be  parabolas  by  examination  of  Equa- 
tion (27b),  which  describes  the  spatial  dependence  of  the  heat  source  selected 
for  these  computations.)  These  computations  were  performed  on  a spatial  mesh 
of  I = J = 31;  the  tick  marks  along  the  boundary  of  the  enclosure  show  the  mesh 
cell  spacing. 

In  Figure  4 isotherms  at  dimensionless  time  t = 10.0  are  shov7n.  By  this  time 
the  flow-field  is  nonlinear,  and  the  temperature  profiles  are  severely  distort- 
ed due  to  buoyancy  effects.  The  temperature  has  increased  and  the  density  has 
decreased  where  heating  has  occurred.  The  heated  fluid  has  become  lighter  than 
its  surrounding  and  begins  to  rise  due  to  buoyancy.  By  continuity,  sur- 
rounding fluid  begins  to  be  drawn  into  the  region  of  the  heat  source,  and  the 
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isotherms,  therefore,  appear  to  be  pinched  off  at  the  bottom  (near  the  center 
of  the  heat  source).  Two  vortices  of  equal  and  opposite  strength  located  on 
the  two  halves  of  the  heat  source  have  developed  to  the  size  where  their  com- 
bined flow  field  is  significant. 

At  time  11.5,  Figure  5,  a buoyant  thermal  has  developed,  giving  the  appear- 
ance of  a mushroom  cloud.  The  two  vortices  mentioned  above  have  begun  to  rise 
with  the  fluid,  being  convected  out  of  the  region  of  primary  heating.  The 
buoyant  thermal  Intensifies  in  strength  as  shown  in  the  next  two  plots  at  times 
12.5  and  13.5,  Figures  6 and  7,  until  the  thermal  hits  the  ceiling,  as  shown 
in  Figure  8,  time  14.5,  and  begins  to  spread.  Inside  the  plume  a distinctly 
periodic  structure  has  begun  to  develop,  as  can  be  seen  vividly  in  Figure  7; 
here,  progressing  up  the  plume  along  its  centerline,  one  finds  a local  low 
first,  then  a periodic  sequence  of  local  highs  and  lows  up  to  about  the  center 
of  the  head  of  the  thermal. 

The  heated  gases  are  seen  to  spread  along  the  ceiling  in  Figure  9 (time  15.5) 
and  fill  the  room  from  the  top  down,  as  shox-m  in  Figures  10,  11,  and  12.  This 
physical  behavior  is  exactly  what  is  observed  in  room-fire  tests  and  in  other 
experimental  observations  of  heating  in  enclosures.  The  S3nnmetry  about  the 
centerline  of  the  room  displayed  in  these  computations  is  some  measure  of  the 
accuracy  with  which  they  were  performed:  the  heat  source  is  placed  S3mmet- 

rically,  but  the  computations  were  performed  as  if  no  symmetry  existed. 

To  assess  the  resolution  of  the  computed  results  shown  in  Figures  3-12,  the 
flow  field  was  computed  again  using  a larger  number  of  mesh  points  I = 63,  vT  = 
64.  Selected  plots  from  this  larger  computation  are  shown  in  Figures  13-18. 
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These  plots  demonstrate  that  the  large-scale  features  determined  by  the  31 
X 31  computation  are  correct  and  agree  with  those  determined  from  the  larger 
computation  to  within  about  ten  percent.  The  results  from  the  larger  computa- 
tion are  characterized  by  smoother  isotherms  and  more  detailed  structures  be- 
cause of  the  greater  resolution.  We  note  that  the  time  required  for  the 
buoyant  thermal  to  reach  the  celling  is  about  ten  percent  longer  in  the 
63  X 64  results  than  in  the  31  x 31  results,  again  apparently  because  of  the 
greater  spatial  (and  temporal)  resolution  of  the  larger  computation. 

Detail  on  a length  scale  of  the  order  of  one  or  two  mesh  cells  must  be  dis- 
regarded because  the  computations  cannot  resolve  such  detail.  On  the  other 
hand,  features  with  a larger  length  scale  can  be  interpreted.  The  spatially 
periodic  behavior  in  the  plume  noted  above  is  a feature  which  requires  some 
discussion.  The  starting  thermal  and  the  plume  induced  by  a heat  source  in  an 
enclosure  are  a result  of  physical  instability  of  the  flow  field.  In  addi- 
tion, in  the  introduction,  we  discussed  the  fact  that  this  fluid  model  was  one 
in  which  viscosity  has  been  ignored,  and,  therefore,  it  could  be  considered  to 
result  from  the  Navier-Stokes  equations  in  the  limit  of  very  large  Grashof 
number  (roughly  Reynolds  number  squared). 

Over  the  last  several  years,  starting  with  the  pioneering  work  of  Brown  and 
Roshko,^®  there  has  been  a reexamination  of  the  meaning  of  turbulence  in  shear 
flows.  There  had  been  a growing  realization  that  turbulence  is  not  satisfac- 
torily described  in  terms  of  velocity  correlations  and  their  corresponding 
spectra  only.  Rather  there  are  distinct  coherent  vortex  structures  in  shear 
flovrs,  and  Brown  and  Roshko  vividly  demonstrated  the  existence  of  these  coher- 
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ent  structures  in  turbulent  shear  flows  using  shadowgraphs  to  visualize  the 
flow  field.  In  particular,  among  many  other  Interesting  features,  Brown  and 
Roshko  found  that  large  scale  coherent  structures  of  the  same  type  existed  in  a 
shear  layer  Independent  of  the  value  of  the  Reynolds  number  provided  only  that 
the  Reynolds  number  is  large  enough  to  have  turbulence.  Subsequent  studies  in 
other  flows,  see  Roshko^^  for  some  of  the  references,  have  shown  that  organized 
structures  exist  in  these  flows  also.  In  addition,  recent  theoretical  studies 
have  shown  that  many  of  the  features  of  the  large  scale  coherent  structures 
observed  in  high  Reynolds  number  flows  can  be  described  by  vortex  structures 
satisfying  Euler's  equations. 

The  spatially  periodic  structures  calculated  in  the  starting  thermal  and  the 
plume  have  been  found  to  be  vortices  of  alternating  sign  produced  in  the  heat- 
ing region  and  convected  out  by  buoyancy.  These  vortices  increase  in  strength 
with  height  above  the  floor  and  occur  in  ant i-S3rmme trie  pairs  with  respect  to 
the  centerline  of  the  room.  We  interpret  these  structures  as  analogous  to  the 
large  scale  coherent  structures  observed  in  turbulent  shear  flows.  In  addi- 
tion, because  these  vortices  are  convected  with  the  buoyant  flov;,  the  spatial 
periodicity  is  translated  into  a temporal  frequency:  at  any  point  within  the 

plume,  each  of  the  dependent  variables  oscillates  with  a frequency  related  to 
the  rate  at  which  vortices  are  generated  and  convected  away.  Experiments,  both 
at  and  elsewhere, have  demonstrated  qualitatively  the  same  feature; 

namely,  buoyant  "puffs",  or  regular  upwellings  followed  by  short  quiescent 
periods,  produced  at  a frequency  determined  by  the  experimental  arrangement. 

The  frequency  of  these  puffs  is  also  found  to  agree  with  the  frequency 
predicted  by  these  calculations. 
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B.  Rectangular  Enclosure  with  Off-Center  Heat  Source 


In  Figure  19  contours  of  constant  temperature  are  shown  at  dimensionless  time 
2.0.  As  in  the  previous  example,  these  contours  mirror  the  contours  of  con- 
stant volumetric  heat  addition  rate  because  the  velocity  field  is  very  small 
and  problem  is  linear  yet  at  this  early  time.  The  heat  source  is  centered 
along  the  floor,  one  quarter  of  the  length  of  the  room  from  the  left  wall. 

At  time  10.0,  Figure  20,  the  Isotherms  have  the  same  form  as  in  the  previous 
example.  Figure  A,  at  time  10.0.  A buoyant  bubble  pinched  at  the  bottom  by 
the  inflow  is  beginning  to  rise  from  the  heat  source.  There  appears  still  to 
be  sjonmetry  about  the  centerline  of  the  heat  source  even  though  the  source  is 
not  symmetrically  placed  within  the  room.  At  time  13.0  a buoyant  thermal  has 
developed,  which  is  asymmetric  with  the  heated  fluid  expanding  more  toward  the 
center  of  the  room  than  toward  the  wall.  Figures  21,  22,  23,  and  24  show  the 
thermal  rising,  growing,  hitting  the  ceiling  and  spreading.  The  heated  gas 
flows  across  the  celling  toward  the  right  in  a gravity  current  while  the  heated 
gas  at  the  left  moves  down  the  left  wall  in  a fashion  similar  to  that  shown  in 
the  first  example  calculation.  As  before,  there  are  large  scale  vortex 
structures  in  the  plume,  but  because  of  the  placement  of  the  heat  source  within 
the  room,  there  is  no  longer  any  symmetry. 
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APPENDIX 


Alternate  Non-dimensional  Variables 


Consider  the  dimensional  system  of  equations  (1)  and  (4)  written  in  the  form: 


3p  3p  1 

+ Ui  + p I j _ 

dt  dX£ 


1 , 

0 - - / Qdv 
V 


0 


= 0 


3ui  3ujL  3p 

p : — + : — + — + [p  - pAy)]  g^i  = 0 


3t 


3x.  3x. 

3 1 


(Al) 


3Ui  y_]^  1 1 

= _[Q--  / QdV] 

3x,  y p V 


Here  p and  uj[  are  the  density  and  velocity  components  as  defined  in  the  text. 
The  quantities  pQ(y)  and  p are  respectively  the  initial  density  stratification 
in  the  vertical  (y)  direction  and  the  difference  between  the  local  pressure  and 
the  hydrostatic  pressure  at  the  height  in  question.  This  pressure  difference, 
which  affects  the  fluid  motion,  is  the  quantity  p defined  in  Equation  (13b)  of 
the  text.  The  heat  source  0 is  prescribed  in  the  form 
'*^0f(t)  1 )2  2 _ 

Q = e c X ^ y (A2) 

^x^y^z  /it 


Note  the  slight  difference  in  notation  from  Equation  (27b)  of  the  text. 

We  now  seek  to  introduce  non-dimensional  variables,  denoted  with  an  asterisk 
(*),  that  are  close  to  those  defined  in  the  text  but  which  reflect  the  strength 
of  the  heat  source  . To  this  end,  we  define  the  following  quantities: 


Xi  = ui  = Uui*(xb*,  t*) 


t 


(H/U)t*.  „ 

> H 

0 


Pap  *(t*) 
0 
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P = Pq(o)  p*(xi(.*,  t*) 


(A3) 


p = Pq(o)  {1  + e p*(X|^*,  t)> 


pQ^y)  = U + 3 pQ*(y*)} 


Here,  is  the  undisturbed  ambient  pressure,  Pq(°)  the  ambient  density  at  the 
floor,  and  H the  height  of  the  enclosure.  The  reference  velocity  U and  the 
dimensionless  density  ratio  3 are  as  yet  undefined.  These  two  quantities  are 
now  determined  by  requiring  that  in  the  Boussinesq  limit,  when  the  density 
ratio  3^0,  all  non-georaetric  parameters  disappear  from  the  problem.  This 
leads  to  the  following  equations  for  U and  3^ 


This  yields  a velocity  scale  which  differs  from  that  employed  in  the  text  by  a 


U^/g  = B 


(AA) 


3U/H  = Qq/H^J^z  ^a* 


/^,  and  a dimensionless  density  ratio  given  by 


2/3 


(A5) 


Finally,  the  equations  of  motion  (Al)  become: 


+ (1  + 3 P*)  D*  = 0 
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(A6) 


Note  that  when  3 is  0(1);  i.e.,  when  there  are  significant  density  variations, 
the  non-dimensionalization  is  for  all  practical  purposes  the  same  as  that  used 
in  the  text.  The  most  significant  feature  of  this  derivation  is  Equation  (A5) 
which  determines  3,  and  hence  the  conditions  under  which  a non-Boussinesq 
model  is  necessary. 
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Figure  Captions 


Figure  1 
Figure  2a 


Figure  2b 
Figure  3 

Figure  4 
Figure  5 
Figure  6 
Figure  7 
Figure  8 
Figure  9 
Figure  10 
Figure  11 
Figure  12 
Figure  13 


Schematic  diagram  of  an  enclosure  with  a two-dimensional  heat 
source,  plume  and  region  of  smoke  and  hot  gases. 

1 

Rectangular  enclosure  in  dimensionless  variables  0 < x < — , 

AR 

0 < y < 1.  The  mesh  upon  which  the  difference  scheme  is  based 
as  shown  schematically  for  (I  = J = 4)  as  a grid  of  solid  lines. 
The  mesh  of  dashed  lines  joins  the  center  points  of  the  basic 
mesh  cells  and  is  the  grid  upon  which  the  pressure  computation 
is  performed. 

A typical  mesh  cell,  with  center  located  at  x = (1  - 1/2)  6x  and 
y = (j  - 1/2)  6y,  Illustrating  where  all  dependent  variables  for 
the  finite  difference  scheme  are  defined. 

Contours  of  constant  temperature  at  dimensionless  time  t = 2.0  in 
a square  enclosure  using  a 31  x 31  mesh.  At  this  early  time  con- 
vection is  unimportant,  and  isotherms  reflect  contours  of  constant 
volumetric  heat  addition. 

Contours  of  constant  temperature  at  dimensionless  time  t = 10.0  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 11.5  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 12.5  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 13.5  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 14.5  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 15.5  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 16.5  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 17.5  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 18.5  in 
a square  enclosure  using  a 31  x 31  mesh. 

Contours  of  constant  temperature  at  dimensionless  time  t = 2.0  in 
a square  enclosure  using  a 63  x 64  mesh. 
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Figure 
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Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

10.0  in 

a square 

enclosure  using  a 63  x i 

64  mesh. 

Figure 

15 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

11.225  in 

a square 

enclosure  using  a 63  x i 

64  mesh. 

Figure 

16 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

13.225  in 

a square 

enclosure  using  a 63  x i 

64  mesh. 

Figure 

17 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

15.225  in 

a square 

enclosure  using  a 63  x 

64  mesh. 

Figure 

18 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

16.187  in 

a square 

enclosure  using  a 63  x 

64  mesh. 

Figure 

19 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

2.0  in  a 

rectangular 

enclosure  using  a 62 

X 31  mesh. 

Figure 

20 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

10.0  in  a 

rectangular 

enclosure  using  a 62 

X 31  mesh. 

Figure 

21 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

13.0  in  a 

rectangular 

enclosure  using  a 62 

X 31  mesh. 

Figure 

22 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

15.45  in 

a rectangular  enclosure  using  a 

62  X 31  mesh. 

Figure 

23 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

17.95  in 

a rectangular  enclosure  using  a 

62  X 31  mesh. 

Figure 

24 

Contours 

of 

constant  temperature 

at  dimensionless 

time 

t 

= 

20.45  in 

a rectangular  enclosure  using  a 62  x 31  mesh. 
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